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ABSTRACT 



In-Situ measurements of the speed of sound in the upper ocean have 
revealed the existence of significant dispersion and large fluctuations 
over the frequency range 25-80 kHz. The near-surface values of c ranged 
from +6 M/sec to -3 M/sec relative to the bubble-free value, with a 
maximum estimated error of 0.5 M/sec. It was possible to identify 
bubbles of "surface" radius centered around 54 microns down to 4.3 meter 
depth as well as a population centered around 124 microns (at 4.3 M) 
found at all depths. The speed fluctuation showed near-Gaussian 
probability density functions except at the dispersion center frequen- 
cies. The standard deviation of the speed varied from 0.27 M/sec for 
58.0 kHz to 0.52 M/sec at 69.6 kHz. The phase remained temporally 
correlated only for 1.46 sec., which was very close to the correlation 
times of the oceanographic variables. The power spectral density of 
the varying phase showed its strongest values at frequencies less than 
0,5 Hz.; this was presumably due to bubbles entrained at orbital frequen- 
cies associated with the surface wave system which had maximum energy 
in the same frequency range. 
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I. ^INTRODUCTION 



The objective of this research was to take in-situ acoustical 
measurements in the surface layer of the ocean in order to determine 
the existence of dispersion and fluctuations of the speed of sound. 

The anticipated source of frequency-dependent speed of sound in sea-water 
are bubbles, of various radii, of numbers dependent on the depth. A 
knowledge of this dispersive behaviour is of importance for sound 
propagation in near-surface ducts where a change in the speed of sound 
with depth would result in. a change in the refraction of the sound ray. 

If bubbles exist near the sea surface a sound ray would bounce more 
often from the surface over a constant horizontal distance as compared 
with the bubble-free case. 

The problem has already been treated in theory as in the analysis 
worked out by Meyer and Skudrzyk (1950). Fox, Curley and Larson have 
examined in the laboratory environment the phase velocity in water 
containing a high density of air bubbles (1954). Glotov et al . (1962) 
studied a population of bubbles generated by wind-agitated breaking 
waves in a laboratory tank. Of the few known data actually taken in 
the ocean those given by Buxcey et al . (1965) and analyzed in the open 
literature by Medwin (1970) will be referenced here. 

Since virtually no in-situ ocean measurements of bubble effects were 
available, this research was proposed and carried out using the U. S. 

Navy Undersea Research and Development Center (NUC) Oceanographic 
Research Tower in Mission Bay/San Diego as a test platform. 

A stable, accurately known, continuous wave acoustical signal was 
received by two spatially separated hydrophones mounted on the sound 



axis and their output was compared in phase by means of a digital phase 
meter. The phase difference between the two hydrophones was recorded 
at each frequency. The frequency was incremented in steps of approxi- 
mately two or four kHz over a range from 25 to 80 kHz and for four 
different water depths. The average results have been interpreted in 
terms of a change in the speed of sound with frequency at a constant 
depth, and with depth at a constant frequency. The statistics of the 
phase change with time have also been calculated. 

At the same time, and at the same location, the local sound amplitude 
modulation due to the spatial and temporal changes in the medium has 
been examined (thesis by LCDR W. J. Smith, December 1971). In addition, 
field measurements were taken to determine the fluctuating temperature, 
salinity, turbulent velocity in x, y, and z-directions , speed of sound 
as given by a velocimeter and wave height in order to allow a better 
judgment of the acoustic-oceanic interaction (theses by Lt. M. Bordy, 

LCDR C. Duchock and Lt. H. Seymour, March 1972), 



II. INSTRUMENTATION 



A steel-pipe frame shaped like a "square C" of height and length 
six ft was used to mount the transducer and the hydrophones as shown 
in Fig. 1. To prevent the frame from acting like a huge tuning fork 
and thereby displacing the hydrophones, the frame was put under tension 
halfway down and at the open end by spring-loaded wires. To minimize 
the reflections from the frame caused by side-lobes of the transducer 
the upper and lower members were covered with acoustic absorbent rubber 
(SoAB) at the specular points. 

Opposite to the open end of the frame and at the center an USRD F27 
unidirectional transducer was mounted. At the sound axis of this trans- 
ducer an Atlantic Research LC-10 hydrophone was placed 78 cm from the 
first one and a total of 173 cm from the F27. This distance between the 
hydrophones was chosen to make the measuring path long enough to measure 
the effect of the bubbles while at the same time limiting the overall 
size of the frame. 

The transmitting part of the instrumentation (Fig. 2) consisted of 
the highly stable General Radio Coherent Decade Frequency Synthesizer 
Type 1162-A generating a sinusoidal voltage of two volts rms and variabl 
in frequency. This signal was then amplified by a Hewlett Packard Power 
467A Amplifier to 14,1 V rms and then impressed across the transducer 
emitting the sound wave. The receiving part of the instrumentation 
(see Fig. 2) consisted of hydrophone #1 with a sensitivity of -125 dB 
re 1 volt/microbar and hydrophone #2 with a sensitivity of -106 dB re 
1 volt/microbar at 40 kHz which picked up the acoustical signal. The 
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Fig. 2. Transmitting and Receiving Instrumentation 




two signals were amplified by 30 dB by an NU5 Pre-Amplifier Model 
2010-030 before being sent through (150) ft of waterproof shielded 
cable. Both signals were then bandpass-filtered by a digitally tuned 
variable Krohn-Hite Model 3322 Filter. The received signal was always 
bandpass-filtered with the center frequency equal to the frequency of 
the transmitted signal through a passband of 1 300 Hz, with -48 dB per 
octave attenuation on either side of the passband. After the Bandpass- 
Filter the signals were fed into a Wiltron Digital Phase Meter Model 355 
which measures the phase angle between the two ac voltages of the same 
frequency and provides a visual display as well as a dc voltage output 
proportional to the measured phase angle at 10 mV per degree. This dc 
voltage was then FM recorded on a Sangamo Model 3500 fourteen track 
magnetic tape recorder simultaneously with the oceanographic data at 
a recording speed of 1 7/8 inch pe.r sec. The dynamic range of the 
recorder utilized was .01 to 2,5 volts, with all the readings falling 
within this range. 



III. OCEANOGRAPHIC COMPONENTS 



It was desired to define the volume of water carrying the acoustical 
signal in terms of the measurable parameters temperature, salinity, 
particle velocity, speed of sound as given by a velocimeter, and surface 
wave height. Therefore ocean sensors were mounted close to the acoustic 
system but not interfering with its performance (see Fig. 3 for geometry). 
Three thermistors (T1 , T2, T3) measured the temporal variation in the 
temperature, and the temperature probe (T4) of the Bissett-Berman 
apparatus gave absolute readings between 14° and 19° C. The Bissett- 
Berman salinometer yielded the temporal variation in the salinity (S) 
for values ranging from 37.5 to 39.5 ppt. The speed of sound c = 
c( T, S, depth) was taken by a Ramsey SVTD probe that provided a temporal 
variation around a DC-value. The particale velocity was obtained with 
an electromagnetic flowmeter reading one horizontal component (u) or (v) 
and the vertical component (w) . A pressure wave height indicator and 
a Baylor gauge gave information on the instantaneous wave height. 
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Spatial Configuration of Acoustic and Ocean Sensors 



IV. EXPERIMENTAL SETUP 



A continuous wave method was used instead of a pulse-echo system 
in order to avoid possible errors related to pulse rise time, definition 
of the leading edge of the pulse, distortion of the pulse when passing 
through the medium because of frequency-dependent attenuation effects 
and the difficulty to determine the phase shift on reflection (see Ref. 
7). This concern was particularly true since we wanted measurements 
over a wide frequency range. 

The method used to measure the speed of sound is actually derived 
from a phase interference approach. The frequency of the signal applied 
to the transducer is adjusted by means of the frequency synthesizer 
until a wanted phase angle is approached at the output of the digital 
phase meter. Since the Sangamo recorder required a signal of at least 
100 mV and at the same time the digital phase meter is most accurate in 
its reading for phase angles close to zero, the frequency was always 
adjusted to get a phase value of approximately 20 degrees. Since the 
digital phase meter compares, the signal at input B with respect to the 
signal at input A and cannot detect integer number of wave lengths 
difference between the two signals, the displayed phase determines the 
fractional wave length by which signal B leads or lags signal A depending 
on the sign. This means that an integer number of wave lengths plus 
this fractional part just fits into the separation distance between the 
two hydrophones. Further increasing the frequency will finally change 
the phase by 360 degrees or step up the integer number of wave lengths 
by one and thus produce the next frequency to be recorded. 



Because bubble populations were expected to be large in the radius 
range 40-150 microns we examined the frequency range from 25-80 kHz. 

The 78. cm separation between the two hydrophones caused the number 
of wave lengths within this spacing to be incremented by one approxi- 
mately every two kHz. In order to exclude possible errors in the results 
due to temporal changes in the medium it was necessary to examine the 
total frequency range in as short a period of time as possible. There- 
fore it was decided to accept 20 minutes as an upper limit for 
stationarity of the medium. The time was also the lower limit to find 
the frequency spectrum of the ocean wave height and turbulent velocity. 

This time allowed us to examine the frequency range in steps of four 
kHz. and to record the time-variable phase difference between the hydro- 
phones for one minute at those frequencies. 

In addition to examining the speed of sound with changing signal 
frequency at a constant depth, readings were taken of the variation of 
the speed for sound with varying depths at a constant frequency. Four 
different depths were chosen in order to observe near-surface to near- 
bottom effects: 4,3, 9.3, 14.3, 7.3 m. These depths were chosen to 

minimize interference from the underwater beam structure of the tower. 

The order was picked to decrease and to check on the possibility of 
longtime variations in the medium. The depths are in meters below mean 
water, where mean water was 6.7 m below the concrete deck of the NUC tower. 

All the readings were taken within 4 1/2 consecutive hours at night 
and early in the morning on October 22, 1971. The obtained data are 
given in Appendix A. 



V. DATA REDUCTION 



Five twenty-minute runs (split into ten two-minute individual 
records) make up the recorded analog data. The approach used was to 
read this information out on a strip-chart recorder, to digitize the 
printed analog data, and to use a computer to carry out a statistical 
and spectral analysis. 

The phase fluctuation and the other oceanographic parameters were 
initially recorded at 1 7/8 inches per second. The tape recorded data 
were played back at 7 1/2 inches per second and first recorded on an 
eight channel Brush Mark 200 strip-chart recorder in order to allow 
visual comparison between the various sensor outputs (Fig. 4). 

The tape recorded data was then played back a second time at 7 1/2 
inches and recorded on a two channel Brush strip-chart recorder. This 
compressed the 20 minute run to five minutes and resulted in an increase 
in frequency by a factor of four. The strip-chart recorder was set at 
25 mm/sec. chart speed and 10 mV/div sensitivity, corresponding to 
1 degree/di v. 

The strip-chart records for each frequency at each depth were digitized 
using the Fleet Numerical Weather Central facility tracing digitizer at 
Point Pinos, California. The digitizer sampled the trace in increments 
of hundredths of an inch in the y-direction for each one one-hundredth 
of an inch advance in the x-direction and recorded this on a nine-track 
magnetic tape. Every record was digitized by first following the zero 
Volt reference line on the strip-chart and then moving the tracer 
vertica-lly at the start of a record from the base-line onto the trace. 



TURBULENT VELOCITY 
X-COMPONENT 




16 




\ 






By doing this it was assured that the digitized data points were 
absolute values referenced to zero Volt. Thus, the scaling factor in 
the y-direction is 

.254 mm/data point = .254 degree/data point 

resulting in 

.254 mm/d.p. X ^ sec/mm = .01016 brush sec/d. p. 

= ,04064 real sec. /data point 

This implies that the sampling rate is 25 Hz and the Nyquist rate is 
then 12.5 Hz. 

The digitized data were read out from the magnetic tape and punched 
on data cards using the Fleet Numerical Weather Central's CDC 5000 
computer (see Computer Program 1). The data cards could then be used 
for analysis on the USNPGS IBM 360 digital computer. 



VI. COMPUTATIONAL METHOD 



Since we were working with a continuous wave method the signal 
received by the first hydrophone can be written as 

y A = a sin (cot) 

and the signal received by the second hydrophone 8 has then to be written 
y B = b sin (ut-kx) 

O JO 

where k = the wave number = 

c c 

f = the frequency of the applied signal 
c = the phase wave velocity 
x = the separation between the hydrophones 

Since we compare the two signals with respect to phase by using the 
digital phase meter we have to find values for co for which 
kx = ip 

where ip is the total phase difference between the two signals due to 
an integral number of wave lengths and a fraction of the wave length. 

We write 

kx = C2Ti)n + 4 > n = an integer 

kx - • ' x = (2 tt)p + <f> 

ax 

(_2tt )n + )> 



c 



fx 




where <}> is in radians 



The final result will be 
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fx 



n + 






360 



where <f> is in degrees 



Therefore the speed of sound may be obtained by using the measured 
frequency at which a phase difference <f> was noted. In addition, x and 
n must be known. The integer n can be found by calculating c from the 
Wilson equation (.K'insler and Frey, Eq. 15.1) using readings for the 
temperature and the salinity. We find 



and then 



VII. BUBBLE THEORY 



The speed of sound through a fluid is given by 



where u> = sound frequency 

k = 2tt/a = propagation constant in bubble-free water 

Following the analysis carried out by Morse and Feshbach, an incident 
acoustic pressure, of wave length very large compared to the bubble 
radius R, having a velocity potential , will produce a total 

field averaged over all possible sizes and configurations of bubbles 
when written in differential equation formulation of 

V 2 < * > + k b 2 (r) O) = 0 



where 
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r = coordinate position of scattering region 

k^ = propagation constant in bubbly water 

R = bubble radius 

u> o = resonant frequency of bubble 

$ = damping constant of bubble 

n(R)dR = number of bubbles per unit volume in radius increment 
R to R + dR 

The quantity — ^ — is the index of refraction in the bubbly region 
and the real part defines the local speed of sound. This term is 
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u(R) dR = fractional air to water volume for the increment between 
R and R + dR 

u(R.)dR. = fractional air to water volume for the increment between 
1 1 R. and R^ + 1/n 



Therefore we can write 
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This is the dispersion relation that predicts a decrease of sound speed 

caused by bubbles of resonance frequency greater than the incident 

frequency and an increase of sound speed for bubbles with lov/er resonant 

* 

frequencies than the impinging sound. 



The resonance frequency of clean air bubbles in water is given by 
(Eckart, 1945) 

f 0 = w 0 / 27T = (1/2 ttR o )£ (3 yP q / p q ) (g/ 

where y = ratio of specific heats of bubble gas (=1.402 for air) 

P Q = ambient pressure at bubble depth 
p Q = density of surrounding water 

The factor g is due to the influence of surface tension, it differs from 
unity only for very small bubbles, for the smallest bubbles of this 
study (R q = 50 microns) g = 1.02. The factor cT^ modifies the value 
of y, i.e,, ycT 1 is the effective ratio of specific heats and lies for 
the smallest bubbles of this study in the range 1,24 < ya~^ 1.40. 

The two effects are to some extent mutually counteracting in their 
effect on the resonance frequency and therefore the term yg/a will 
be assumed to equal 1.33. 

The damping factor as a function of bubble resonance frequency 

as given by Devin will be approximated by 

0,322 

= ,027 Iw 0 ./27i x 10 3J 
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VIII. EXPERIMENTAL RESULTS 



Five different freqency sweeps have been carried out, runs 6 and 7 
at a depth of 4,3 m and the remaining three, runs 8, 9, and 10, at 
9,3, 14.3 and 7.3 m depth, respectively. 

Run 6 covered a total frequency range from 25 - 80 kHz in steps of 
approximately two kHz. to be able to detect the detailed variations 
in the speed of sound; for the remaining runs readings were taken every 
four kHz, in order to examine the frequencies of interest in as short 
a period of time as possible. Run 7 started on hour after the end of 
run 5, but for the five runs reported here all the readings were taken 
within 4 1/2 consecutive hours between 0350 and 0820 on October 22, 1971. 

Single temperature and salinity readings taken at each depth and 
a single bathythermograph record taken at 8 a.m. (Fig. 5) allow one 
to calculate the speed of sound for bubble-free water; the Wilson 
equation is used in the following form: 



c = 1449 + 4 ,6t - 0.055t 2 +0.0003t 3 + (1 .39-,0l2t) (S-35) + 0.017d 



where c = speed of sound in meters per second 
t = temperature in degrees centigrade 
S = salinity in parts per thousand 
d = depth in meters 

The bubble-free values obtained for 4.3, 9.3 and 14.3 m are used as 
references for the dispersion curves and are included in Fig. 6. 

The calculated values for the speed of sound based on the obtained 
data are examined and analyzed with respect to a frequency — and/or depth 
dependence. If there is any change in c with frequency then this change 
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is attributed in general to either molecular relaxation or bubbles 
since for the sea-water medium no other frequency -dependent mechanisms 
are known at the frequencies of this study, But the magnitude of the 
molecular relaxation effect can be shown to be negligible (a frequency 
shift from 20 to 70 kHz. results in a Ac~0.3 m/sec.) as compared to 
the bubble effects . 

A, FREQUENCY BEHAVIOUR 

All the calculated values of c as a function of frequency for the 
different runs have been plotted on Fig. 6 for comparison. They should 
be interpreted by making use of the reference value for the non-bubbly 
medium which, because of the temperature change, shifts by as much as 
4.2 m/sec. in going from 4.3 m to a depth of 14.3 m. 

The first point of interest is the dispersive variation of c at a 
single depth. The first run (6) was carried out at a depth of 4.3 m 
and examined frequencies from 30.7 up to 79.2 kHz during a 40 minute 
period of time. Most significant is the strong increase in the speed 
of sound by as much as 6 m/sec. compared to the non-bubbly value for 
c of 1510.47 m/sec. and caused by the peak at 34.7 kHz followed by a 
smaller negative change of 3 m/sec at 40.6 kHz. For frequencies from 
47 to 60 kHz the speed of sound shows nearly no dispersive behaviour. 
The higher frequencies are characterized by another peak in the c-curve 
roughly one fourth the strength of the peak at 34.7 kHz. 

Run seven is interesting since it covered the same frequency range 
and the same depth as run six but for the major part now in steps of 
four kHz rather than 2 kHz. Since this allows a comparison between the 
two runs, a reminder is necessary that there is a time lag of one hour 
between the end of run six and the start of run seven. This time delay 



may account for the slight shift in the dispersion curve downward and 
to the right. The shape matches the previous one pretty well, with two 
peaks clearly discernible. 

In run eight the depth is increased by five meters as compared with 
the previous run. The main feature is a resonance pattern with a weaker 
maximum at 36.8 kHz rather than 34.7 kHz at shallower depth. 

Run nine, at the greatest depth (14.3 m) , differs in the pattern 
for the speed of sound over the frequency range by showing an increase 
in c over a wide frequency band with a maximum at approximately 41 kHz 
and no second peak. The obvious shift in the level is caused by a 
decrease of the non-dispersive value of c from 1510.47 at 4.3 m to 
1506.27 m/sec at 14.3 m. 

Run ten was made back at a depth of 7.3 m. The two readings at 
57,7 and 61.59 kHz are believed to be non-representative since a boat 
approaching the tower injected a large number of bubbles into the volume 
under consideration. Excluding these data points, the dispersion curve 
shows a main peak around 35.5 kHz. 

Summarizing the dispersion data it can be stated that: 

1) The change of the speed of sound with frequency at sea was observed 
to be significantly stronger than expected. 

2) The change in the speed of sound does not follow the usual dispersive 
pattern for a single dominant bubble regime. 

3) The speed of sound curve at 4,3 m has two peaks. The frequency of 

the lower peak changes from 34.7 kHz at 4.3 m to 41 kHz at 14.3 m, which 
is a shift toward higher frequencies as expected for surface-generated 
bubbles carried downwards "isothermal ly" by turbulence. However, the 
shift is smaller than expected: isothermally convected identifiable, 



non-diffusing, surface-initiated bubbles which resonate at 35 kHz at 
4.3 m would resonate at 54 kHz and not at the observed 41 kHz when 
they reach 14.3 m. There are two explanations that can be proposed 
for this discrepancy: 

a) the presence of the higher frequency peak in the near-surface 
measurement has shifted the low frequency peak to higher frequencies 
than it would have had by itself. The absence of the second peak in 
the deeper measurement (due to more rapid disappearance of smaller 
bubbles) permits the low frequency peak to more accurately identify 
the resonance frequency of these bubbles. We herein assume that near 
the surface we encountered a bubble population made up of larger surface- 
generated bubbles and smaller bubbles introduced primarily by continental 
aerosols that drop into the water. To check on the dispersion curve that 
results from a specific bubble distribution, a bubble population was 
created peaking around 60 and 120 microns as given in Fig, 7. Subse- 
quently the computer program (3) was written calculating the speed of 
sound dispersion based on the equations derived in Section VII and using 
the above mentioned volume concentration of bubbles. The result is a 
dispersion curve given in Fig. 8 which is in good agreement with our 
data obtained at 4.3 m. The observed dispersion curves for the other 
depths support this hypothesis, since only part of the larger bubbles 
are carried down, decreasing in number with increasing depth whereas the 
smaller bubbles rapidly disappear due to gas diffusion out of the 

bubbl es . 

b) the larger bubbles originate in the volume or at the bottom, rather 
than at the surface. As they rise they lose gas by diffusion, so that 
when they reach the surface they are not as large as predicted by the 
non-diffusion assumption and are thereby at higher frequencies than expected. 
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Speed of Sound Dispersion for Bubble Population of Fig. 



B. DEPTH BEHAVIOUR 

Next the variation in c with depth for a constant frequency will 
be considered, On the assumption that the bubble densities were not 
a function of time during the 2 1/2 hours of runs 7—10, all these data 
have been rearranged to display the variation in the speed of sound 
with depth at various constant frequencies as given in Table I. These 
values have been used to plot the speed of sound versus depth in Fig, 9. 
These straight line segments reflect the major variations in c with 
frequency as described in the last section. 

In order to examine the variation of the speed of sound versus 
depth, it is necessary to subtract out first all possible temperature, 
salinity, and pressure-based changes with depth. To do this, use is 
made of the following coefficients for variation with 

temperature || = + 5 ft/ (sec ) C°F) 

salinity || = + 4 ft/ (sec) (ppt) 

depth || = + 0.017 ft/(sec) (ft) 

The salinity changed by -.1 ppt over 10 m and therefore accounts for a 
change of c with S of -.4 meter per second. It is necessary to add to 
this a change of +.17 m/sec due to the depth effect such that the total 
effect of salinity and depth on c will be -.23 m/sec when going from 
4.3 to 14.3 m; this small effect will be neglected. 

The temperature-caused changes in the speed of sound have to be 
taken into consideration. The BT curve sketched in Fig. 5 is approxi- 
mated by straight line elements and is used with the Wilson equation 
to calculate values for the speed of sound at different depths for the 



TABLE I. Speed of Sound in m/sec 



^Depth^Ori) 
Freq. (kHz) 


4.3 


4.3 


7.3 


9.3 


14.3 


28.97 








1505.73 


1503.28 


30.70 


1511 .99 


1511 .37 


1507.67 






32.77 


1514.18 


1513.88 




1509.89 




34.76 


1516.46 


1516.94 


1512.51 






36.79 


1512.90 


1516.52 




1511 .55 


1509.2 


38.73 


1509.48 


1511.82 


1510.93 






40.61 


1507.30 


1510.22 




1509.89 




42.50 


1508.13 


1509.30 


1509.53 






44.47 




1509.39 




1510.17 


1509.89 


46.40 


1510.29 




1509.62 






48.36 


1511 .31 


1510.05 




1509.84 




50.27 


1511 .36 




1508.78 






52*20 


1511 *03 


1509.79 




1509.36 




54.1 2 


1510.74 




1509.20 






56,03 


1510.92 


1509.99 




1510.66 


1506.56 


59*88 


1510.31 


1509.19 




1510.64 


1506.78 


63,67 


1509.11 


1507.65 




1508.95 


1505.43 


65.62 


1508.11 




1507.86 






67.54 


1509.14 


1507.55 




1507.92 




71.42 


1501 .32 


1509.26 




1508.42 


1505.01 



Speed of Sound versus Depth 
for various frequencies (in kHz) 
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the non-bubble case. After doing this it is possible to compare the 
measured values of c with the values for bubble free sea water. The 
plot of the temperature-caused changes has to be referenced to a value 
of c that holds for the bubble-free medium and the previous calculated 
value of 1506,27 m/sec at 14.3 m was used. The non-bubble theoretical 
value of c is subtracted from the measured value and the difference is 
then called ac^, due to bubbles, which is plotted on the same graph. 

Since the shift in the frequency peaks causes the pattern of the measured 
c versus depth plot to be markedly different for the various frequencies, 
five different frequencies have been picked that are felt to be repre- 
sentative for the total frequency region: 36.79, 44.47, 56.03, 63.67, 

71,42 kHz. 

This representation of c versus depth at a constant frequency can 
only be used to examine whether one encounters, at a specific depth, 
resonant bubbles such that for sound frequencies smaller than the bubble 
resonance frequency one experiences a decrease of sound speed and for 
impinging sound of higher frequencies than the resonance frequency an 
increase in c. It should be clear that this plot ignores depth effects 
on an identifiable surface-generated bubble carried to different depths. 

Figure 10 presents this information for a frequency of 36.79 kHz. 
Since there is a strong positive Ac^ of 5.25 m/sec at 4.3 m this means 
that we encounter strong populations of bubbles with lower resonance 
frequencies than the impinging sound, in other words, there are many 
bubbles that resonate below 36.79 kHz. Figure 11 represents the situ- 
ation at 44.47 kHz and shows that we sense bubbles of resonance frequency 
greater than the incident frequency at 4.3 m depth and of resonance 
frequency smaller than the incident frequency at 14,3 m. The picture 
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changes when going to Fig. 12. At a frequency of 56.03 kHz there is 
no significant dispersive behaviour at all. For frequency 63.67 kHz, 
Fig. 13, at all depths we sense bubbles of resonance frequencies greater 
than the frequency of the acoustical signal. Since the difference 
between the measured value and the non-bubble theoretical value is 
much larger at 4.3 than at 14.3 m we are either very close in frequency 
to resonant bubbles at 4.3 m or the population is considerably stronger 
than at 14,3 m. Figure 14, 71.43 kHz, finally shows a smaller negative 
ACg with only minor differences between the various depths. 

Summarizing, we can make the following statements related to the 
different depths: 

4.3 m : There is a large population of bubbles having a resonance 

frequency close to and somewhat smal ler than 36.79 kHz. At 44.47 and at 

56.03 kHz we sense and at 63.67 kHz we detect another population of 
bubbles with an ay greater than 63.67 kHz. 

7.3m : Bubbles with a resonance frequency less than 36.79 kHz are also 

detected at this depth. At 44.47 kHz we are far enough away from this 
region and encounter no dispersive characteristic. 

9.3 m : There is a population of bubbles resonating below 36.79 kHz. 

14.3 m : There are bubbles of resonance frequencies below 44.47 kHz. 

At 56.03 kHz this region of resonance frequencies has been passed and 
no dispersive behaviour is shown. 

Conclusions : 

At 4,3, 7.3, and 9.3 m we observed bubbles of resonance frequency below 
36.79 kHz and above 71.42 kHz and at 14.3 m we found a broad region of 
bubbles resonating below 44.47 kHz. The sound speed dispersion is 
strongest at 4,3 m and is reduced in amplitude when going to a greater 



depth. The conclusions derived from the variation in the speed of sound 
with depth are consistent with the observations made for the variation 
in the speed of sound with frequency. 



IX, STATISTICAL AND SPECTRAL ANALYSIS 



The signal received at hydrophone one does not have a constant phase 
difference with respect to that at hydrophone two but is a random signal. 
This is the equivalent of a noisy signal where the noise in this case is 
due to the integrated effects of time-varying speed of sound due to 
varying temperature, salinity and bubble patches along the sound path. 

It is assumed that the random signal exhibits statistical stationarity 
which is partially described in terms of a probability density function 
(P.D.F.) and the statistical moments. In addition to a statistical 
description it is valuable to use the time series to determine the 
temporal correlation function and the pov/er spectral densities of the 
time-varying phase (Ref. '2). 

The variable phase signal that we received at a specific frequency 
and for a period of approximately 90 sec is a continuously changing 
random variable within some finite interval. But since this recorded 
information is digitized, the random variable can assume only a finite 
number of distinct values so that we must treat it as a discrete random 
variable. In practice the recording technique used resulted in magni- 
tude increments of .254° spaced in time by .04064 sec. For one frequency 
record we then get n data points. If A is one of the possible data points, 
then n^ is the number of times A occurred and the ratio n^/n is the 
relative frequency of occurrence of the event A for that particular 
record. This approach allowed the use of the sampled data points to 
form a histogram that can be interpreted as a distribution function: 

F(x) = PU^x) 



from which a density function can be obtained: 



p(x) = £ P(x.) 5(x - x.) 
i 

Examining these histogramsr it turns out that the main pattern of the 
envelope follows a Gaussian distribution. This can be expected since 
in the experimental measurement many irregular and fluctuating causes 
like weak temperature and salinity inhomogeneities caused the phase 
reading to have small random variations which produce the Gaussian 
P.D.F. about the undisturbed value. It seems that the bell-shaped 
curve experiences stronger changes around 35 and 65 kHz. This different 
behaviour is pointed out by a series of historgrams (Figs. 15-19) for 
the latter frequency range and for 4.3 m depth. At 59.9 kHz the general 
shape is narrow Gaussian but at 63. 7 kHz the distribution is nearly 
evenly spread over a wider center range that sharply drops at the 
skirts. At 57.6 kHz the envelope changes again and achieves at 69.6 
kHz a wide Gaussian pattern. At 77.3 kHz the pattern is again that of 
a narrow Gaussian noise. Since this change in the density function 
is a particular phenomenon for frequencies around 35 and 65 kHz and is 
not found for frequencies in between, it is postulated that time varying 
resonant bubble populations predominantly around frequencies 35 and 65 
kHz caused the variation with frequency. These frequencies are essen- 
tially the same as observed previously in the peaks in the analysis of 
speed as a function of frequency. 

The statistical average, or mean value, 3T of the random variable <f> 
has been used to calculate the mean speed of sound which is the expected 

value for c. By calculating the variance 

2 T =72 ~Z - 2 

<? x = (x - x) = x - x 
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Fig. 17. Histogram of Phase Angle 
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and taking the square root of the variance to obtain the standard 
deviation, c , one has a measure of the spread of the observed values. 

A large standard deviation implies a large spread of likely values for 
any given observation, and vice versa. The standard deviation was 
calculated for all frequency readings from 30.7 to 79.2 kHz and is 
plotted as a function of frequency in Fig. 20. 

Since there was a time lag of two minutes between frequency runs, 
a is for this case not only a function of frequency but also of time, 
which is increasing with frequency. To be able to judge which fraction 
of the variation in a is a time effect and which is the frequency 

A 

effect a continuous 20 min record of phase modulation at 60 kHz was 
split into ten two-minute records so that the standard deviation could 
be examined as a function of time alone (Fig. 21). Comparing the two 
graphs there is a mean value of 1.28° and a total variation of i .32° 
for the time record; there is a much stronger excursion for the frequency 
plot which shows a mean value of 1.92° at 60 kHz with a deviation of 
+ 1.92° at 40.6 and of + 2.56° at 69.6 kHz. One can rule out temperature 
and salinity patches as the causing agency for the effect in Fig. 20 
since those would not show a frequency dependent spread of the measured 
values. It seems that the presence of bubbles in the medium causes this 
strong change in the standard deviation as a function of frequency. The 
frequencies of the peak values of a are approximately the same as the 
frequencies already identified as centers of resonant bubbles. 

Next the autocorrelation function R(t) will be considered since it 
provides information about <^(t) and also relates to the frequency-domain 
description of the random signal. R(t) is defined by' 

R(t) = < x(t)x(t+x) > / < x(t)x(t+x)> 

x=0 
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Fig. 20. Standard Deviation as a Function of Frequency. 
Time increasing with frequency 




Fig, 21. Standard Deviation as a Function of 
Time Only at Frequency 60 kHz 



It should be noted that the autocorrelation function is a measure of 
both temporal variation and statistical dependence. The plots of the 
normalized autocorrelation at 30.7 and 46.4 kHz are considered to be 
representative and are given in Figs. 22 and 23. The graph for 30.7 
kHz shows a variation superimposed on the almost ideal damped cosine 
pattern - obtained at 46.4 kHz. This variation is probably due to bubble 
effects. Recall that the frequency 30.7 kHz was observed as a predo- 
minant resonance frequency in the speed of sound discussion of section 
VIII. The time it takes for the autocorrelation to drop to 1/e of its 
value is plotted versus frequency on Fig. 24, No specific regularity 
can be found. On the average this "correlation time" is approx. 1.8 sec. 

Examining the autocorrelation function we arrive at the following 
results: The phase modulation remains correlated (1/e) for 1.46 sec, 

the salinity for 3.96 sec and the speed of sound as measured by the 
velocimeter for 2.28 sec (Figs. 25-27). The temperature records have 
not yet been analyzed (thesis of LCDR Duchock to be issued March 1972). 
Since for the velocimeter (operating at 3 MHz) c = c (T, S) it appears 
that the temporal variation in c is due principally to salinity patches 
(and most probably temperature patches). On the other hand, our phase 
measurement is a function of not only S and T but also of the bubble 
distribution. The still smaller correlation time of 1.46 sec is evidence 
of additional decorrelation due to bubbles. 

Turning to the frequency domain, the power spectral density of the 
phase variation gives the distribution of average power of the fluctu- 
ations in frequency. To find the PSD the Mean Lag Product Method by 
B1 ackman/Tukey was used. The main feature of this method is to calculate 
not only the autocovariance function but in addition an apparent 
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Fig. 22. Autocorrelation at 30.7 kHz 
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Fig. 23. Autocorrelation at 46.4 kHz 
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Fig. 24. Decorrelation Time as a Function of Frequency 
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Fig. 25 



Autocorrelation for Phase Modulation at £§ kHz 
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Fig. 26. Autocorrelation for Speed of Sound Fluctuations 
by Velocimeter 
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Fig. 27. Autocorrelation for Salinity Fluctuations 



autocovariance by using a Parzen-lag-window. By taking the Fourier 
transform of this function one gets a smoothed power spectrum. The 
result of this computation using computer program (2) is a plot of phase 
angle squared per Hz versus frequency (The frequency resolution is .05 
Hz. For better resolution in the amplitude, 10 log-^ of the ordinate 
is taken such that the dimension is in dB). Several conclusions can 
be derived from the PSD-plot for the 4.3 m depth. For nearly all 
frequencies there is a distinct peak around .48 Hz probably due to the 
orbital effects caused by the observed predominant ocean wave component. 
From 0.5 to 2.5 Hz the values for the PSD follow an exponentially 
decaying curve. This pattern is interrupted by peaks mainly at .8, 

1,3, 1,55, 1.95, and 2.35 Hz that may show up at two or at most three 
consecutive readings (see peak at 2.3 Hz in Figs. 28-31), or occur at 
one frequency and the next time two records later. Since every recording 
took two minutes this then suggests that the agency forcing these changes 
in the frequency spectrum is of a nature that can result in variations 
that show up only over a time period of four minutes or repeat every 
four minutes. 

In general temperature, salinity, and/or bubble patches carried by 
turbulence and crossing the acoustical path could be the driving force 
for these variations. Since the sound signal integrates all these 
effects it is necessary to examine the temperature and salinity records. 
For the 20 min run and the depth of 4.3 m the power spectral density 
of the phase modulation with a frequency resolution of .0045 Hz, of the 
Ramsay velocimeter reading, and of the Bissett-Berman salinometer, all 
follow approximately the same exponentially decaying pattern, Figs. 

31-33 (these graphs courtesy of LCDR Duchock and LT Seymour). The 
principal difference is that the PSD of the phase shows more power 
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Fig. 28, Power Spectral Density, 46.4 kHz 
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Fig. 29. Power Spectral Density, 48.4 kHz 
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Ftg. 30. Power Spectral Density, 50.3 kHz 
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Fig. 31. Power Spectral Density for 
Phase Modulation at 60 kHz 
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Fig. 32. Power Spectral Density for Speed of Sound 
Fluctuations by Velocimeter 
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Ftg. 33. Power Spectral Density for 
Salinity Fluctuations 



concentrated below 0.5 Hz where orbital velocities are strong. All 



three plots of 
the range from 



PSD show good agreement even in the fine structure 1 
0.5 - 1.0 Hz. 



X. ERROR ANALYSIS 



The speed of sound as derived before is given by 



c = ^ = _I2L 



n ± <f / 360 



In order to define the accuracy with which c is given we use the 
following approach: 

log c = log f + log x - log (ri ± <j>/360) 

dc _ df + dx^ _ d(n ± <j>/360 ) 

c f x n + ^/360 

then the maximum error is given by 



This error will be calculated for the highest (80 kHz) and the lowest 
(30 kHz) frequencies. The coherent decade frequency synthesizer 
synthesizes the frequency of the output signal by proper combination 
of a number of frequencies, each derived on a proportional basis from 
a single internal master frequency crystal. It was essential in this 
application because of the requirement for a stable, sine-wave signal 
source capable of ultraprecise frequency adjustments through a wide 
frequency range. Since the output is always coherent in frequency with 
that of the quartz-crystal oscillator source it has the same percentage 
accuracy as the source which is essentially the smallest digit in the 

_7 

specific frequency or 1.6 x 10 at 60 kHz. This then implies that 



Ac 

c 




A(n ± 4>/ 360 
n- ± <f> / 360 



Af 

T 



k = 3.3 x 10 

30x1 O' 3 



-7 



at 30 kHz and 



* 








Af 

f 



at 80 kHz 



10 



-2 



80x10" 



= 1 .25 x 10 



-7 



Next the ax/x term will be examined. The absolute hydrophone distance 

O 

could only be measured to one part in 1.5 x 10 (0.5 mm) which is 

equivalent to 



AX 

X 



5x10 

.78 



-4 



= 6.4 x 10 



-4 



Assuming the Af/f and the A ip/\p terms to be zero this error in measuring 
the separation between the hydrophones results in an absolute error for 
the speed of sound of 1 m/sec (a level shift of 1 m/sec in the non- 
bubbly value of c), but does not affect the dispersion curve, that is 
the pattern of the speed of sound versus frequency. Therefore this 
error is no serious problem and will be ignored in the analysis since 
we are not interested in the absolute accuracy but rather the dispersion 
in the speed of sound . 

To reduce possible errors due to changes in the geometrical config- 
uration, especially in the separation distance between the hydrophones, 
it was necessary to use as a mounting a wire put under 100 lb of tension. 
The remaining possible vibration due to forces at sea was estimated to 
be 0.1 mm and results in 



AX 

X 




1.28 x 10~ 4 



The accuracy of the phase was the limiting factor. Since basically 
a phase interference approach was used to measure the speed of sound 
any relative phase changes could not be allowed to occur in the two arms 
of the system. It is a question whether one can assume the two parallel 



parts within the system (hydrophone-preamplifier-filter-input channel 
phase meter) to be identical with respect to their influence on the 
signal. Actually, complete identity can never be obtained over a wide 
frequency range. Therefore special attention was given to the phase 
characteristic in selecting the hydrophones. The LC-5 resp. 

LC-10 were selected because their resonance is far above the frequency 
range of interest. Comparing the dispersion curve for 4.3 m with that 
at 14.3 m also points out that there is no phase variation in the 
pattern that shows up consistently at all depths and may therefore be 
due to a different phase behaviour of the two hydrophones. In addition 
the zero control of the phase meter has been used to compensate at a 
single frequency for unequal phase in the two receiving parts when the 
two input signals were known to be exactly in phase by using the 
electrical signal from the output of the frequency synthesizer. 

The accuracy of the analog output of the phase meter is ±0.1°±0.3% 
of the phase angle measured for all frequencies or ±0.16° for the 20° 
phase angle which was used normally in the experiment. This causes a 
fractional error 



A(n ± 360 ) = .16 

n ± <J>/360 ' (360 x 16) + 20 



= .272 x 10~ 4 at 30 kHz 



.16 



(360 x 41) + 20 



.108 x 10~ 4 at 80 kHz 



There are of course accuracy limitations in measuring the phase of 
distorted wave forms. Because the phase meter makes its phase measure- 
ments by measuring the fraction of a cycle between the negative-going 
zero axis crossing of the two test signals harmonic or scattered signal 



can displace the zero crossing in either a plus or a minus manner 
depending on the phase relationship to the test signal. 

The scattered signal was measured to be down by 37.4 dB re 1 V/y bar 
as compared with the main signal which for the worst case could result 
in an error of .77°, A scattered signal a adding to the main signal b 
would result in the 1 argest error in e when the condition depicted 
below occurs 



such that 

for small e 
then 

and 

A possible source of error is also the fact that at some frequencies 
the response of one hydrophone is significantly stronger than the other 
which means the possibility of harmonic content in the signal at 
multiples or submultiples of these frequencies. The limit error for 
a harmonic of magnitude, e, is 6 = — 57° where a is the magnitude of 

cl 

the fundamental (manual for digital input phase meter model 355). 

Since the second harmonic was measured to be -46.34 dB re 1 V/ybar 
as compared with the fundamental the harmonic content could at most 
result in an error of 0.28°. 

Assuming the worst condition for the phase reading or an error of 
.16° at the phase meter, a maximum additional error of .77° due to the 




0 = = .0135 in radians 

2 ,77° 



scattered signal, and also a maximum error of .28° due to harmonic 
content in the signal the total phase error could be 



A(n ± 4>/360) 
n + <j)/360 



.16 + .77 + .28 
(360 x 16) + 20 



2.08 x 10~ 4 at 30 kHz 



1 .207 

(360 x 41 ) + 20 



.82 x 10‘ 4 at 80 kHz 



Summarizing it can be stated 

— = 3.3 x 10~ 7 + 1.28 x 10~ 4 + 2.08 x 10" 4 = 3.36 x 10" 4 at 30 kHz 

c 

= 1.25>xl0~ 7 + 1.28 x 10“ 4 + .82 x 10" 4 = 2.1 x 10~ 4 at 80 kHz 



Therefore the maximum error in c is mainly determined by the percentage 
error in x and in and is given by 

Ac = (1500) (3.36 x 10“ 4 ) = .5 m/sec at 30 kHz and 

= (1500) ( 2.1 x 10' 4 ) = .31 m/sec at 80 kHz 

Comparing the value for Ac at the lov/er and at the upper limit points 
out that in general one is able to take more accurate measurements at 
the higher frequencies. 



XI. SUMMARY AND CONCLUSIONS 



A fixed source-hydrophone system separated by 1.8 m was used to 
study a) the in-si tu speed of sound as a function of frequency and 
depth and b) the temporal characteristics of the phase fluctuations 
of a continuous wave acoustical signal over the frequency range 25 - 80 
kHz. The experiment was carried out during Sea State 1 conditions at 
the NUC Oceanographic Research Tower in water of depth 56 ft one mile 
off shore at Mission Bay/San Diego in October 1971. The near surface 
region is of particular interest because bubbles and turbulence can 
have significant acoustical effects in addition to those caused by 
temperature and salinity inhomogeneities. 

Four different depths, 4.3, 7.3, 9.3, and 14.3 m, were examined. 

The near surface values of c ranged from + 6 m/sec to -3 m/sec relative 
to the value for bubble-free sea water, with a maximum error of 0.5 m/sec. 
The speed of sound was found to peak at approximately 35 and 71 kHz at 
4.3 m depth. The lower peak frequency increased as the depth is in- 
creased but not at the rate which would be expected for bubbles that 
are wind-initiated at the surface and carried to greater depths without 
losing their identity. A possible explanation is that the bubbles 
originated within the volume or at the bottom and lost gas by diffusion 
as they rose to the surface. Another possibility is that presence of 
more high frequency bubbles at the surface affects the apparent 
resonance frequency of the larger bubbles as determined by sound dis- 
persion. The gradual and smooth drop-off in the resonant behaviour with 
increasing depth demonstrates that bubbles of this average size are 
somewhat continuously distributed in the medium down to 14.3 m. The 



second distinct peak centered at 71 kHz appears to be only a character- 
istic of the shallowest run and does not repeat at the other depths. 

It is postulated that these small bubbles are introduced at the surface 
either by randomly breaking waves or by continental aerosols that drop 
into the water and carry tiny air bubbles with them (Medwin, 1970). 

The surface initiated bubbles disappear quite fast due to gas diffusion 
out of the bubble when transported downwards in the medium. We conclude 
from the foregoing that we encountered bubbles of radii centered around 
124 and 54 microns and resonating at 30 and 69 kHz respectively. The 
broad dispersive characteristic around 30 kHz implies many bubbles 
larger than 60 microns whose principal source was at the surface and 
which were entrained by convective currents carrying them downwards. 

Statistical analysis shows that for frequencies other than at the 
peak bubble populations the phase fluctuations can be approximated by 
a Gaussian probability density function. However around the peak 
variations in sound speed the PDF loses the smooth Gaussian appearance 
and shows a wide, somewhat irregular spread. The temporal variation 
of the standard deviation of the sound phase was small (± .32°). How- 
ever, the standard deviation of the phase fluctuation changed 
considerably as a function of frequency increasing to 1.92° at 40.6 
kHz and to 2.56° at 70 kHz, the two sound dispersion centers of this 
experiment. 

Typically, the phase was decorrelated 0/ e ) after 1.46 sec. This 
compares with preliminary oceanographic information which shows corre- 
lation times of 3.96 sec for the salinity and 2.28 sec for the speed 
of sound as measured by the velocimeter. These values show internal 
consistency. 



The power spectral density of the varying phase shows its strongest 
values at frequencies less than 0.5 Hz; this is presumably due to 
bubbles entrained at orbital frequencies associated with the surface 
wave system which has maximum energy in the same frequency range. 

Due to time considerations it was not possible to further examine 
the phase fluctuations in terms of the functional dependence on 
temperature, salinity, turbulence, and wave height. The analysis of 
the oceanographic parameters will be carried on in the theses by 
Lt. Bordy, LCDR Duchock, and Lt. Seymour (March 1972) and the emphasis 
should then be placed on the interrelation of the parameters as given 
by the cross-correlations or cross spectra for instance. 



XII. RECOMMENDATIONS 



It is recommended that research in this area be continued, especially 
by making an in-situ determination of the bubble concentrations at the 
same time as sound speed measurements are taken in order to determine 
the variation of bubble population with location, v/eather conditions, 
day or night, and local effects. It is also suggested to examine a 
smaller frequency range in more detail rather than to skim through a 
wide frequency range in larger steps and thereby losing valuable 
information. Frequencies below 30 kHz seem to be most inviting because 
of relevance to the Navy and promising because the largest changes in 
the speed of sound occur here, just below the frequency of the most 
prominent bubble population. 



APPENDIX A. PRESENTATION OF DATA 



Run 6, 4.3 m 



freq 

[kHz] 


real t 
[sec] 


mean 

[Dig.] 


standard 

dev. 

[Dig.] 


mean 

[radians] 


n 


f X 

c=~ x=,783 


30.72 


180 


-.84 


.06 


-.07467 


15.0987 


1511 .99 


32.70 


101 


-.83 


.08 


-.07378 


16.9096 


1514.18 


34.71 


90 


-.69 


.06 


-.06133 


17.9220 


1516.46 


36.82 


80 


+ .445 


.11 


+.03956 


19.0562 


1512.90 


38.70 


30 


+ .65 


.08 


+.05778 


20.0745 


1509.48 


40.60 


38 


+ .83 


.12 


+.07378 


21.0905 


1507.30 


42.42 


44 


+ .08 


.08 


+.00711 


22.0238 


1508.13 


46.40 


193 


+ .44 


.06 


+.03911 


24.0558 


1510.29 


48.40 


107 


+ .665 


.06 


+.05911 


25.0758 


1511 .31 


50.30 


102 


+ ,48 


.06 


+.04267 


26.0593 


1511.36 


52.22 


96 


+ .485 


.03 


+.04311 


27.0598 


1511 .03 


54.14 


102 


+ .49 


.06 


+.04356 


28.0602 


1510.74 


56.09 


99 


+ .57 


,06 


+.05067 


29.0673 


1510.92 


58.00 


98 


+ .49 


.06 


+ .04356 


30.0602 


1510.77 


59.92 


90 


+ .54 


.06 


+.04800 


31 .0647 


1510.31 


61 .82 


89 


+ .53 


.06 


+.04711 


32.0638 


1509.65 


63.72 


89 


+ .50 


,07 


+,04444 


33.0611 


1509.11 


65.64 


103 


+ .51 


,11 


+.04533 


34.0620 


1508.90 


67.59 


121 


+ .58 


.11 


+.05156 


35.0682 


1509.14 


69.60 


104 


+ .66 


.14 


+.05867 


36.0753 


1510.64 


71 .56 


86 


+ .65 


,10 


+,05778 


37.0745 


1511 .32 


73,49 


90 


+ .70 


,11 


+,06222 


38.0789 


1511.14 


75.38 


101 


+ .59 


.09 


+,05244 


39.0691 


1510.72 


77.28 


110 


+ . 56 


.07 


+.04978 


40.0665 


1510.25 


79.20 


78 


+ .57 


.07 


+.05067 


41 .0673 


1510.05 



Run 7, 4.3 m 



freq 

[kHz] 


real t 
[sec] 


mean 

[Dig.] 


standard 

dev, 

[Dig.] 


30.69 


73.6 


-.94 


.06 


32.69 


112.8 


-.85 


.04 


34.69 


113.6 


-.87 


.05 


36.60 


88.0 


-.97 


,04 


38.76 


98.4 


+ .65 


.09 


40.61 


93.6 


+ .43 


.07 


40.61 


93.6 


+ .44 


,07 


42,53 


104.8 


+ .53 


.05 


44.45 


93.6 


+ .47 


.04 


48.34 


88.8 


+ .55 


.05 


52.19 


92.8 


+ .56 


.06 


56,04 


104.8 


+ .48 


.06 


59.86 


78.4 


+ .45 


.09 


63.66 


99.2 


+ .51 


.09 


67.50 


99.2 


+ .47 


,07 


71,44 


73.6 


+ .52 


.09 


75.24 


92.0 


+ .50 


.08 


79.10 


105,6 


+ .51 


.07 



mean n c= — x=.783 

[radians] n 



-.0836 


15.8997 


1511.37 


-.0756 


16.9077 


1513.88 


-.0773 


17.9060 


1516.94 


-.0862 


18.8971 


1516.52 


+.0578 


20.0745 


1511 .82 


+.0382 


21.0549 


1510.22 


+.0391 


21 .0558 


1510.16 


+.0471 


22.0638 


1509.30 


+.0418 


23.0585 


1509.39 


+.0489 


25.0656 


1510.05 


+.0498 


27.2667 


1509.79 


+.0427 


29.0594 


1509.99 


+.0400 


31 .0567 


1509.19 


+.0453 


33.0620 


1507.65 


+.0418 


35.0585 


1507.55 


+.0462 


37.0629 


1509.26 


+.0444 


39.0611 


1508.22 


+.0453 


41.0620 


1508.34 



Run 8, 9.3 m 



freq 

[kHz] 


real t 
[sec] 


mean 

[Dig.] 


standard 

dev 

[Dig.] 


mean 

[radians] 


n 


c=— x=.783 
n 


25.17 


81.6 


+ .51 


.03 


+.0453 


13.062 


1508.81 


28.98 


79.2 


+ .60 


.03 


+.0533 


15.070 


1505.73 


32.91 


84.8 


+ .56 


.05 


+.0498 


17.0665 


1509.89 


36.79 


96,0 


+ .46 


.03 


+.0409 


19.0576 


1511 .55 


40.62 


91 .2 


+ .54 


,03 


+.0480 


21 .0647 


1509.89 


44.49 


93.6 


+ .57 


.06 


+.0507 


23.0674 


1510.17 


48.34 


102.4 


+ .59 


.07 


+,0524 


25.0691 


1509.84 


52,18 


95.2 


+ .59 


.04 


+.0524 


27.0691 


1509.36 


56.07 


103.2 


+ .51 


..10 


+,0453 


29.0620 


1510.66 


59.94 


88.8 


+ .58 


.08 


+.0516 


31 .0683 


1510.64 


63.71 


101 .6 


+ .48 


.08 


+.0427 


33.0594 


1508.95 


67,56 


108.0 


+ .49 


.07 


+.0436 


35.0603 


1507.92 


71 .39 


54.4 


+ .46 


.05 


+.0409 


37.0576 


1508.42 



Run 9, 


14.3 m 






freq 

[kHz] 


real t 
[sec] 


mean 

[Dig.] 


standa rd 
dev, 
[Dig.] 


28.95 


42.4 


+ .70 


.05 


36.76 


28.0 


+ .62 


,08 


44.48 


32.0 


+ .56 


.04 


55.93 


36.8 


+ .58 


.06 


59,78 


36.8 


+ .54 


.05 


63.57 


25.6 


+ .53 


,06 


71,28 


78.4 


+ .76 


.09 



mean 

[radians] 


n 


c=^ x= , 783 


+.0622 


15.0789 


1503.28 


+.0551 


19.0718 


1509.20 


+.0498 


23.0665 


1509.89 


+ .0516 


29.0683 


1506.56 


+.0480 


31.0647 


1506.78 


+.0471 


33.0638 


1505.43 


+.0676 


37.0843 


1505,01 



Run 10, 


7.3 m 






freq 

[kHz] 


real t 
[sec] 


mean 

[Dig.] 


standard 

dev, 

[Dig.] 


27,05 


81.6 


.59 


.04 


30,96 


71.2 


.70 


.05 


34,89 


96,8 


.51 


.04 


38,72 


89,6 


.55 


.04 


42,55 


87.2 


.61 


.04 


46.40 


89.6 


.56 


,05 


50.23 


81 .6 


.57 


,09 


54,09 


86,4 


.52 


,05 


57.70 


84,8 


*-1 ,08 


.07 


61.59 


80,0 


-.79 


.06 


65,59 


78,4 


.48 


..13 



mean 

[radians] 


n 


c=^ x=.783 


.0524 


14.0691 


1505.44 


.0622 


16.0789 


1507.67 


.0453 


18.0620 


1512.51 


.0489 


20.0656 


1510.93 


.0542 


22.0709 


1509.53 


.0498 


24.0665 


1509.62 


.0507 


26.0674 


1508.78 


.0462 


28.0629 


1509.20 


-.0960 


29.8873 


1511.65 


-.0702 


31.9131 


1511.13 


.0427 


34.0594 


1507.86 



oooooono ooooo 



PROGRAM 1 

PROGRAM CONVERT ( INPUT, OUTPUT , PUNCH) 

AVGX IS THE AVERAGE DISTANCE DISTANCE THE RECORDING 
DRUM ADVANCES PER ONE HOUR REAL TIME, I.E. THE DATA 
SAMPLING RATE, 

DIMENSION U(6000)»V(6000)»N(80)»IBUFF( 5000),NK(80) 

I K=0 

READ 97, L 

97 FORMAT (12) 

READ 98, AVGX 

98 FORMAT ( F10.0) 

READ 96 , 1 F ILE 

96 FORMAT (12) 

DO 116 J J= 1 , L 
DELT=36 . 0/ AVGX 
PRINT 99 , DELT 

99 FORMAT ( 1H0,2X,25H SAMPLING INTERVAL EQUALS, 1X,E15.7, 
11X,7HSEC0NDS) 

PRINT 700, JJ 

700 FORMAT (1H0,9HJJ EQUALS, 2X,I2> 

DO 100 1=1,5500 
U( I ) =0 . 0 

100 V(I)=0.0 

101 DO 202 1=1,5000 
202 I BUFF ( I ) =0 

COUNTX=0 . 0 
C0UNTY=0 . 0 
NUM=5000 
M= 1 
K8=0 

CALL L IOF( 5LRBCD1, I BUF F, NUM, NPAR , NEOF ) 

IF(NEOF) 602,200 

200 K=-7 
KF=0 

201 K=K+8 
KA=K+8 
KC=0 

DO 102 KB=K, KA 
k r — k r + i 

NK(KC)=IBUFF(KB) 

102 IF (NK(KC) .EQ.O) KF=1 

DECODE (80,103,NK) (N(I),I=1,80) 

103 FORMAT ( 80R 1 ) 

DO 104 19=1,80 

IF(N( 19) .EQ.50B) GO TO 106 
IF(N( 19) .E0.55B) GOTO 107 
I F ( N ( 19) .EQ.349) GO TO 108 

SYMBOL / (50) REPRESENTS AN INCREMENT TRAVEL IN THE 
MINUS X OR Y DIRECTION BY THE DIGITIZER. 

SYMBOL 0 , ( 55B ) , REPRESENTS A ZERO INCREMENT TRAVEL 
IN THE X OR Y DIRECTION BY THE DIGITIZER 
SYMBOL 1,(345), REPRESENTS AND INCREMENT TRAVEL IN 
THE POSITIVE X OR Y DIRECTION BY THE DIGITIZER 
SYMBOL *,(473), IS A FLAG INSERTED IN IN THE RECORD 
BY THE PERSON DIGITIZING BY USE OF THE IRG 
GO TO 104 

106 RX=-0 .01 
K8=K8+1 
GO TO 109 

107 RX=0 . 0 
K8=K8+ 1 
GO TO 109 

108 RX=0 . 0 1 
K8=K8+ 1 

109 K3=K8/2 
K3=2*K3 

IF(K3.EQ.K8) GO TO 111 
COUNTX=COUNTX+RX 
I F ( COUNT X. NF .0 . 0 ) GO TO 110 
GO TO 104 

110 U ( M ) =COUNTX 
GO TO 104 



oo 



111 COUNTY=COUNT Y+RX 

I F { COUNTX . NE .0 . 0 ) GO TO 112 
GO TO 104 

112 V ( M ) =C0UNT Y 
C0UNTX=0 .0 
M=M+ 1 

I F ( M .GT . 6000 ) GO TO 600 

104 CONTINUE 
IF(KF.EQ.l) GO TO 113 
GO TO 201 

113 MAX=M— 7 

105 PRINT 205, M, 19 
205 FORMAT ( 1H0 » 21 1 0 ) 

TOTAL TIME OF THE RECORD EQUALS NUMBER OF DATA 
POINTS TIMES THE SAMPLING INTERVAL, DELT. 
TIME=(M*DELT)/3600 .0 
PRINT 115, TIME 

115 FORMAT !1H,10X,28H TOTAL TIME OF RECORD EQUALS, 2X, 
1E15.7,2X,7H HOURS. ) 

PRINT 215, ( V( I ) , 1=1 ,M) 

PRINT 215, (U( I) , 1=1, M) 

215 FORMAT ( 1H ,10X,14F7.2> 

PUNCH 213, ( V( I ) , 1=1, M) 

213 FORMAT ( 14F5 .2 ) 

PUNCH 2100 

2100 FORMAT ( 20H******************-* ) 

116 CONTINUE 
STOP 

600 PRINT 601 

601 FORMAT! 1 HO, 20X ,37H***** U AND V SPACE INADEQUATE ***,) 

602 IK=IK+1 

I F ( IK.GT.IFILE ) 603 , 200 

603 STOP 
END 



no o oonnonnononoonnnnnno 



PROGRAM 2 

PROGRAM FOR THE SPECTRAL ANALYSIS OF APERIODIC RECORDS 
DATA MUST BE SAMPLED AT EQUAL INTERVALS OF DT 
CALCULATES EITHER THE AUTO SPECTRA OR THE CROSSSPECTRA 

IF NFL AG1=0 , AUTOSPECTRA IF NFLAG1 = 1 COSPECTRA 

IF NFL AGP = 1 CALL PRESS, TOTE 

IF NFLAGC = 1 CALL URMS 

NTS= NUMBER OF TIME SAMPLES 
MLAG * NUMBER OF TIME LAGS 
DT = TIME INCREMENT 
DHZ = FREQUECY SPACING 
FBHZ = LOWEST FREQUENCY (HZ) 

FEHZ = HIGHEST FREQUENCY OF INTEREST (HZ) 

TM =DT*MLAG 
DHZ = 1/(2*A*TM) 

FEHZ = B/(2*TM)= B/ ( 2*MLAG*DT ) 

FEHZ = A*B*DHZ 

A*B = TOTAL NO. FREQUENCY SPACINGS OF ENERGY SPECTRUM 
A AND B ARE INTEGER CONSTANTS 
TRANS — INSTRUMENT AND REMARKS 

REAL*8 ITITLE( 12) , ITITELt 12) , ITITIL(12) 

DIMENSION C 1(6000),F2( 6000 ) 

DIMENSION PHI (600) ,TAU(600) ,SPHI (600) ,APHI (600) , 

1PHN( 600) 

DIMENSION FREQ (700) , CYCL ( 700 ) , SPEC ( 700 ), S SP ( 700 ) 
DIMENSION CSPEC(700) ,QSPEC(700 ) ,SS2(700) 

DIMENSION RESPF( 700) , PER (700), SPE2( 700) , PHASE (700) , 
1C0HER(700) 

DIMENSION DATE (2), HOUR (2), TRANS (10) 

REAL*4 LABEL/4H CO / , L ASL I /4HQUAD/ , MONT , LABI L/4H PHI/, 
1LABLE/4HCR0S/ 

READ! 5,200, END=888) ( ITITLE( I ) , 1=1 ,12 ) 
READ(5,200,END=888) ( ITITEL ( I ) ,1=1,12) 

200 FORMAT ( 6A8 ) 

READ ( 5 , 802 ) DATE , HOUR, TRANS 

802 FORMAT ( 2 A4 , 2 A4 , 1 0 A4 ) 

READ (5,96) NFL AG 1,NFLAGP, NFLAGC 
READ(5, 99) NTS , MLAG , DT , DHZ , F BHZ , FEHZ 
READ (5, 801) CALX 1,CALX2,H,X2,AZZ,S 
WRI TE ( 6,811) 

811 FORMAT ( 1H1 ) 

WRITE( 6 , 803 ) D AT E , HOUR , TRANS 

WRITE(6, 98) NT S, ML AG, DT, DHZ, FBHZ, FEHZ , CALX1 , CALX2 , H, X2 
96 FORMAT (316) 

98 FORMAT ( 7H NTS = ,I6//,8H MLAG = ,I6//,6H DT = ,F9.5//, 
16H DHZ= , F 1 0 

2.8//,* FBHZ = , ,F10.8//,» FEHZ = *,F12.8//,* CALX1 = • 
3,F10.8// , 

4* CALX2 = •,F10.8//,» H = »,F10.3//,' X2 = ',F10.3//) 

99 FORMAT (2I6,5F10.8) 

801 FORMAT (6F10.5) 

803 FORMAT ( 40H SOUND PHASE STUDY OF SAN DIEGO BAY , 

15X,7H DATE , 2 A4 , 1 0 H HOUR , 2 A4// , 5X, 1 0A4/// ) 

BAND WIDTH FREQUENCIES OF TOTAL ENERGY FLUX- CMIN,CMAX 
CMIN = 0.0 
CMAX = 0.2 

COMPUTING POWER SPECTRUM 
PI=3. 14159265 
FB = FBHZ*2.0*PI 
FE = FEHZ*2.0*PI 
DF = DHZ *2.0*PI 
FMIN = 2.0*PI*CMIN 
FMAX = 2.0*PI*CMAX 
900 FORMAT ( 14F5. 2 ) 

READ ( 5 , 900 ) ( F 1 ( I ), 1=1, NTS) 

XMAX = Fl( 1 ) 

DO 23 1=2, NTS 

IF(XMAX.GE.F1( I ) ) GO TO 23 
XMAX = FI ( I ) 

23 CONTINUE 

DO 24 1=1, NTS 



FI C I ) = ( F 1 ( I)-XMAX) 

24 CONTINUE 

WRITE! 6, 901) (FI ( I), 1=1, NTS) 

901 FORMAT ( 14F7. 2 ) 

CALL TREND(F1 ,NTS,DT,CALX1) 

21 DO 10 M= 1 , MLAG 
SUM=0.0 
NMAX=NT S-M+l 
DO 8 1=1 ,NMAX 
NN=M+I-1 

8 SUM=SUM + F1( I )*F1( NN) 

XNMAX=NMAX 
XX=M- 1 

TAU(M)=XX*DT 
PHI ( M) =SUM/XNMAX 
PHN(M) = PHI ( M ) / PH I ( 1 ) 

10 CONTINUE 

CALL PARZ ( MLAG , PHI ) 

798 FORMAT ( ’ PHI(O) = *,F10.4//) 

NFREQ=(FE-FB)/DF+0.1 
DO 14 N=1 , NFREQ 
XN=N 

FREQ (N ) = ( XN-1 . 0 ) *DF+FB 

14 CYCL(N) = FREQ(N)/ (2.0*PI) 

PER(l) = 0.0 

DO 15 N = 2, NFREQ 

15 PER ( N ) = 1.0/CYCL(N) 

MLAGM1=ML AG-1 
XMLAG=ML AG 

46 DO 50 N= It NFREQ 

SUM=0.5*(PHI(1 )+PHI(MLAG)*COS(FREQ(N)*TAU(MLAG) ) ) 
C1=C0S(FREQ(N) *DT) 

S1=SIN(FREQ(N)*DT) 

CC=1.0 

SC=0.0 

DO 49 M=2, MLAGM1 

CT=CC*C1-SC*S1 

ST=SC*C1+CC*S1 

CC=CT 

SC=ST 

49 SUM=SUM+PHI ( M ) *CC 

50 SPEC ( N ) = SUM*2. O/XMLAG 

I F ( NFL AGP ) 741,741,742 
742 CONTINUE 

CALL PRESS (FREQ, SP EC, NFREQ, H , X2 , DF, FB , FMAX) 

741 CONTINUE 

MD=3.14159265/(DF*DT*XMLAG)+0.1 
DC = DF/ ( 2 . 0 *P I ) 

SUM = 0.0 
DO 650 N=l, NFREQ 

650 SUM = SUM + SPEC(N) 

WRITE (6,651) SUM 

651 FORMAT ( 25H VARIANCE OF SP ECTRUM/ / , 3 X , 

116H VARIANCE = ,F10.5,10H M2 //) 

WRITE! 6, 106) ( TAU(M) ,PHN(.M) ,M = 1 ,MLAG) 

INUM=FFHZ/DHZ 

DO 9988 1=1, INUM 

SPEC ( I )=10*AL0G10(SPEC(I ) ) 

9988 CONTINUE 

C CALL DRAW TO PLOT SPECTRUM OF SAN DIEGO 

C AS A FUNCTION OF FREQUENCY 

CALL DR AW ( 200 , T AU , PHN , 0 , 0 , LABEL , I T ITEL , 8 . 0 , 0. 0 , 2, 0 , 2 

1 .0. 8.8. 1, LAST) 

CALL DR A W( I NUM , C YCL , SP EC , 0 , 0 , LAB EL , I T ITL E , 0 .0 , 0. 0 , 0 , 0, 

10. 0. 8. 8.1, LAST) 

105 FORMAT <10X,24H ENERGY SPECTRAL DENS ITY// / , 5X , 

1 1 OH FREQUENCY, 5X,13H 

2 RAW SPECTRUM, 3X, 18H SMOOTHED SPECTRUM//) 

106 FORMAT ( 5X , F 10 . 5 , 6X , F 1 0. 5 ) 

WRITE(6, 105) 

WRITE( 6 , 106 ) (CYCL( M),SPEC(M) , M=l, NFREQ) 

888 RETURN 



END 



SUBROUTINE PRESS (FREQ, SPEC , NFR EQ , H, X2, OF , FB , FMAX ) 

C SUBROUTINE TO CONVERT PRESSURE SPECTRUM TO ENERGY SPECTRUM 
DIMENSION FREQ(NFREQ) , SPEC (NFR EQ ) 

NMAX = ( FMAX— F BJ/DF+1.0 
DO 11 1=1, NMAX 

C CALCULATE LINEAR WAVE LENGTH BY NEWTONS METHOD 
IF(FREQ( I)-0. 00001) 8,8,7 

7 XKHO = FREQ( I)*FREQ(I)*H/9.80 
IF ( XKHO-6. 3 ) 5,1,1 

1 XKH = XKHO 
GO TO 9 

5 XKH = SQRT(XKHO) 

3 SH = SINH(XKH) 

CH = COS H( XKH ) 

EPS = XKHO-XKH*SH/CH 
SLOPE = -XKH/CH**2-SH/CH 
DXKH = -EPS/SLOPE 
IF(ABS(DXKH/XKH) —0 .0001) 9,9,4 

4 XKH = XKH+DXKH 
GO TO 3 

8 RESPF = 1.00 
GO TO 11 

9 XK = XKH/H 

RESPF = C0SH(XK*H)/(C0SH(XK*(H-X2))) 

I F ( RES PF— 1 0 .0 ) 11,11,12 

11 SPEC(I) = RESPF*RESPF*SPEC(I ) 

12 RETURN 
END 



SUBROUTINE TREND ( FX , NTS , DT , C ALXX ) 

DIMENSION FX ( NTS ) 

C CALI BRAT I ONG RECORD 
DO 104 1=1, NTS 
104 FX(I) = FX(I)*CALXX 
C COMPUTING THE LINEAR TREND 
FNTS = NTS 
SUMF = 0.0 
DO 101 1=1, NTS 

101 SUMF = SUMF + FX ( I ) 

SUMF 1 = 0.0 

DO 102 1=1, NTS 
XI = I 

102 SUMF1 = SUMF1 + XI*FX(I) 

XNM1 = NTS-1 

XNP1 = NTS+1 

XM = (1.0/DT)*(12.0*SUMF1/(FNTS*XNM1*XNP1)-6.0*SUMF/ 
1 (XNM1*FNTS) ) 

B = SUMF/FNTS-XM*XNPl*DT/2.0 
FMEAN = SUMF /FNTS 
WRITE ( 6 , 9 ) FMEAN, XM,B 

9 FORMAT ( 3X , 8 H MEAN = ,F10.5,3X,9H SLOPE = ,F10.5,3X, 
1 13H I NTERC EPT= ,F10.5//) 

DO 103 1=1, NTS 
XI = I 

103 FX ( I ) = FX ( I ) - (B+XM*XI*DT) 

RETURN 

END 



SUBROUTINE SMO ( MD, XI , X2 , NFREQ ) 
DIMENSION XI (MD) ,X2(MD) 

DO 1 N= 1 , MD 
NA=N+MD 
NN=NFREQ-N+1 
NB=NN-MD 

X2(N) = 0.25*<X1(1)+X1(NA))+0.5*X1(N) 
1 X2(NN)=0.5*(X1 (NN)+X1(NB) ) 

3 MB=MD+1 



ME=NN-1 

5 DO 2 N=MB» ME 
NA=N+MD 
NB=N-MD 

2 X2(N)=0.25*(X1(NA)+X1(NB) )+0.5*Xl(N) 
RETURN 
END 



SUBROUTINE PARZ ( MLAG, PHI ) 

C PARZ SUBROUTINE PARZEN FILTERS AUTO-CORRELLATION FUNCT 

DIMENSION PHI ( ML AG ) 

XMLAG = MLAG 

MLAGH = XMLAG/2. O-O.l 

MLAGH1 = MLAGH + 1 

DO 31 M=l, MLAGH 

MM = M-l 

R = MM 

RM = R/XMLAG 

UM = 1.0-6.0*RM*RM*(1.0-RM) 

PHI ( M ) = PHI ( M )*UM 

31 CONTINUE 

DO 32 M = MLAGH1 f MLAG 
MM = M-l 
R = MM 

RM = R/XMLAG 
RM1 = ( 1.0-RM) 

UM = 2.0*RM1*RM1*RM1 
PHI ( M ) = PHIIM )*UM 

32 CONTINUE 
RETURN 
END 



C PROGRAM 3 CALCULATES DISPERSION IN C 

// EXEC FORTCLGP, REGION. G0=100K 
//FORT . SYS IN DD * 

REALMS T ITLE ( 1 2 ) / ' BOX 209 ',11*' '/ 

REAL LABEL/' •/ 

DIMENSION A (10) ,B( 1 0 ) , MR ( 1 0 ) , W ( 500 ) , CC ( 500 ) 

DATA PI/3.1415926/ 

C SET MAX DIMENSION TO WORK WITH. 

NOUT = 500 

C READ PZ,RHZ,CZ. 

READ (5,100) PZ,RHZ,CZ 
100 FORMAT ( 3 FI 0 • 5 ) 

PZ=PZ*1Q0000.0 

C READ RADIUS (MICRONS) AND FREQ. 

READ (5, 200) I S R , I DR , I ER , IS W, I DW, I EW 
200 FORMAT (6110) 

C READ NO. OF SEGMENTS AND THEIR VALUEX. • .MAX OF 10 SEGMENTS 
READ ( 5 , 300 ) NS 
300 FORMAT (I 2) 

READ (5, 400) ( A( I ) ,B ( I ) , MR( I ) , 1 = 1 ,NS) 

400 FORMAT ( 2F1 0 . 5 , 1 10) 

C PRINT INPUT. 

WR I TE ( 6 , 500 ) 

500 FORMAT ( • 1INPUT : • ) 

WRIT E( 6 » 510 ) PZ,RHZ,CZ 

510 FORMAT ( * PO = ' , El 6 . 8, T28 , ' RHO = • , E16. 8 ,T54, 'CO = ', 
1E16.8) 

WRITE (6,520) ISR,IDR,IER 

520 FORMAT ( ' STARTING RADIUS = ' »I6»T28,' INCREMENT = ',13, 
1T54, 'ENDING RADIUS = ',16) 

WRITE (6,530) ISW,IDW,IEK 

530 FORMAT ( ' STARTING W = ', 16, T28,' INCREMENT = ',I3,T54, 

1 'ENDING W = ', 16 ) 

WRITE( 6, 540) 

540 FORMAT ( • SEGMENT • ,TX , 'A* , 15X, 'B» ,9X, ' MAX. RADIUS') 
WRITE (6, 550) ( I , A< I ) , B ( I ) , MR ( I ) , I=i ,NS ) 

5 50 FORMAT ( 5X, I2,3X,E16.8»2X,E16.8 ,4X_» 16) 

C CHECK FOR INPUT ERROR. MR SHOULD INCREASE. 

NJ=NS- 1 

IF (NJ .EQ. 0) GO TO 700 
DO 600 J=1 ,NJ 

IF ( MR ( J ) .LT. MR ( J+l ) ) GO TO 600 
WRITE (6,560) 

560 FORMAT (• SEGMENT RADII ARE OUT OF SEQUENCE.') 

STOP 

600 CONTINUE 

C NOW CALCULATE SOME CONSTANTS TO SAVE TIME. 

700 IR=( IER-ISR) /I DR+1 
IW=( IEW-ISW)/IDW+1 
WC = 1 ,975*SQRT( PZ/P.HZ ) * 1000000 . 0 
P I C=0 .027/ ( < P 1^2 00 0. 0)**0. 322) 

C CLEAR OUTPUT BUFFER 
DO 800 1=1, NOUT 
W( I ) =0 . 

800 CC(I)=0. 

C INITIALIZE STARTING "VALUES. 

C WW=WORKING OMEGA. 

C WZI=CALCULATED OMEGA ZERO SUB I. 

C R W=WGRK I NG RADIUS. 

I TEMP W= I SW 
DO 2000 K=1 , I W 
WW=I TEMP W 
W( K)=ITEMPW 
IF (K.NE.l ) GO TO 915 
WRITE16,900) I TEMP W 
900 FORMAT ( • 1W = • , 16) 

WRITE ( 6 , 91 0 ) 

910 FORMAT ( • INDEX* , 6X, • RADI US ', 1 IX , 'WZ' , 13X ,' U ( R) D ( R )• ) 
915 SUM=0 . 0 

ITEMPR= ISR 

C THIS IS THE SUMMATION LOOP. 

DO 1000 L= 1 , IR 



WZI=WC/ ITEMPR 
DI=PIC*(WZI**0.322) 

WRAT=WZI/( WW#2 .0*P I *1000.0) 

WRAT2=WRAT*WRAT 
WRATM=WRAT 2-1.0 

BATH=WR ATM/ ( (WRATM*WRATM)+(DI*D I* W RAT 2 ) ) 

C SELECT U(R)D(R). 

DO 920 J = 1 » NS 

IF (ITEMPR .GT.MR(J)) GO TO 920 
URDR=1.0E-09*( A( J)*ITEMPR+B( J) ) 

BUBBLE=URDR/ ( ITEMPR*ITEMPR*1 .0E-12) 

GO TO 927 
920 CONTINUE 

C IF FALLS THROUGH HERE, WE HAVE GONE TOO FAR. 

925 WR I TE ( 6 , 926 ) 

926 FORMAT!* ERROR ON SEGMENT RADIUS.*) 

STOP 

C WRITE RESULT FOR EACH I. 

927 IF (K.NE.l) G3 TO 940 
WRITE(6,930) L , I TE MPR , WZ I , UROR 

9 30 FORMAT ( 2X , 1 6 , 6X , 16 , 6X , E 16 . 8 , 2X , E16 . 8! 

940 ITEMPR=ITEMPR+IDR 

SUM=SUM+( RUBBL E*BATH ) 

1000 CONTINUE 

C NOW CALCULATE SOUND VELOCITY. 

CNEW=1 . 0/ ( ( 1 .O/CZ) +( (3.0*CZ*SUM)/(8,0*PI*PI*WW*WW* 
11.0E06) ) ) 

CC (K)=CNEW-CZ 

WRITE (6,1010 ITEMPW,SUM,CNEW,CC(K) 

1010 FORMAT ( * W = *,I6,4X,'SUM = • , E 16 . 8,2X , • C = *,E16.8, 
12X, 'DIFF = ' ,E16.8 ) 

C PREPARE FOR NEXT ROUND OF W. 

ITEMPW= I TEMPW+I DW 
2000 CONTINUE 

C ALL DONE.. NOW PLOT W AND C. 

WRITE! 6,2100 ) 

2100 FORMAT ( • 1 RESULT: W VS. ( C— CO ) • ) 

WRITE(6,2200)(W(I) ,CC( I) ,I=1,IW) 

2200 FORMAT! 1 OX , E20.10, 10X, E20. 10) 

WRITE( 6,2300) 

2300 FORMAT Cl') 

CALL PLOTP <W,CC,IW,0) 

STOP 

END 
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